1 



Variational wave functions for frustrated 
magnetic models 

Federico Bccca, 1 Luca Capriotti, 2 Alberto Parola, 3 and Sandro Sorclla 1 

1 CNR-INFM-Democritos National Simulation Centre and International School 
for Advanced Studies (SISSA), Via Beirut 2-4, 1-34014 Trieste, Italy 

2 Global Modelling and Analytics Group, Investment Banking Division, Credit 
Suisse Group, Eleven Madison Avenue, NY-10010-3629 New York, United States 

3 Dipartimento di Fisica e Matematica, Universita dell'Insubria, Via Valleggio 11, 
1-22100 Como, Italy 

Variational wave functions containing electronic pairing and suppressed charge 
fluctuations (i.e., projected BCS states) have been proposed as the paradigm 
for disordered magnetic systems (including spin liquids). Here we discuss the 
general properties of these states in one and two dimensions, and show that 
different quantum phases may be described with high accuracy by the same 
class of variational wave functions, including dimerized and magnetically or- 
dered states. In particular, phases with magnetic order may be obtained from 
a straightforward generalization containing both antifcrromagnetic and super- 
conducting order parameters, as well as suitable spin Jastrow correlations. In 
summary, projected wave functions represent an extremely flexible tool for 
understanding the physics of low-dimensional magnetic systems. 



1.1 Introduction 

The variational approach is a widely used tool to investigate the low-energy 
properties of quantum systems with several active degrees of freedom, includ- 
ing electrons and ions. The basic idea is to construct fully quantum many-body 
states by a physically motivated ansatz. The resulting wave function should be 
simple enough to allow efficient calculations even for large sizes. Most of the 
variational calculations are traditionally based upon mean-field approxima- 
tions, where the many-body wave function is constructed by using independent 
single-particle states. In this respect, even the BCS theory of superconductiv- 
ity belongs in this category [T]. Although these mean-field approaches have 
been instrumental in understanding and describing weakly correlated systems, 
they have proved inadequate whenever the electron-electron interaction domi- 
nates the kinetic energy. The generalization of variational states in this regime 
is not straightforward, and represents an open problem in the modern theory 
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of Condensed Matter. Probably the most celebrated case is the wave function 
proposed by Laughlin to describe the fractional quantum Hall effect as an 
incompressible quantum fluid with fractional excitations [2J. One important 
example in which electron correlations prevent the use of simple, mean-field 
approaches is provided by the so-called resonating valence-bond (RVB) state. 
This intriguing phase, which was conjectured many years ago by Fazekas and 
Anderson [3], has no magnetic order, no broken lattice symmetries, and re- 
mains disordered even at zero temperature. It is now commonly accepted that 
these spin-liquid states may be stabilized in quantum antiferromagnets with 
competing (frustrating) interactions [4]. 

Here we present one possible approach to the definition of accurate varia- 
tional wave functions which take into account quite readily both strong elec- 
tron correlations and the frustrated nature of the lattice. The price to pay 
when considering these effects is that calculations cannot be performed an- 
alytically, and more sophisticated numerical methods, such as the quantum 
Monte Carlo technique, are required. 

Let us begin by considering a simple, frustrated spin model, in which 
the combined effects of a small spin value, reduced dimensionality, and the 
presence of competing interactions could lead to non-magnetic phases. We 
consider what is known as the J1 — J2 frustrated Heisenberg model on a chain 
or a square lattice, 

h = Ji s r ■ s «' + J 2 s * • s «'> ( L1 ) 

n.n. n.n.n. 

where J\ and J2 are the (positive) nearest-neighbor (n.n.) and next-nearest- 
neighbor (n.n.n.) couplings, and Sr = (S^, S^, Sff) are 5 = 1/2 operators; 
periodic boundary conditions are assumed. Besides the purely theoretical in- 
terest, this model is also known to describe the relevant antiferromagnetic in- 
teractions in a variety of quasi-one-dimensional [5] and quasi-two-dimensional 
systems 

In one dimension, the phase diagram of the J1 — J2 model has been well es- 
tablished by analytical studies and by Density Matrix Renormalization Group 
(DMRG) calculations [8]. For small values of the ratio J2/ Ji, the system is in a 
Luttingcr spin-fluid phase with a gapless spectrum, no broken symmetry, and 
power- law spin correlations. By increasing the value of the second- neighbor 
coupling, a gapped phase is stabilized 08]. The value of the critical point has 
been determined with high accuracy as (J 2 /Ji) c = 0.241167 ± 0.000005 0. 
The gapped ground state is two-fold degenerate and spontaneously dimerized, 
and at J2/J1 = 0.5 is expressed by the exact Majumdar- Ghosh wave func- 
tion |10[lll|. Interestingly, for J2/J1 > 0.5, incommensurate but short-range 
spin correlations have been found, whereas the dimer-dimer correlations are 
always commensurate [8]. 

By contrast, the phase diagram of the two dimensional J1 — J2 model is 
the subject of much debate. For J2/J1 <C 0.5, an antiferromagnetic Neel or- 
der with magnetic wave vector Q = (7r,7r) is expected. In the opposite limit, 
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Jil J\ ^ 0.5, the ground state is a collinear antifcrromagnctic phase where 
the spins are aligned ferromagnetically in one direction and antifcrromagnct- 
ically in the other [Q = (ir,0) or Q — (0,7r)]. The nature of the ground state 
in the regime of strong frustration, i.e., for J2/J1 ~ 0.5, remains an open 
problem, and there is no general consensus on its characterization. Since the 
work of Chandra and Doucot [12] , it has been suggested that a non-magnetic 
phase should be present around J 2 1 ' Ji — 0.5. Unfortunately, exact diagonaliza- 
tion calculations are limited to small clusters which cannot provide definitive 
answers to this very delicate problem [T3J[T3J[T5] . By using series-expansion 
methods [16l[T7l[T8l[T9] and field-theoretical approaches [20], it has been ar- 
gued that a valence-bond solid, with columnar dimer order and spontaneous 
symmetry-breaking, could be stabilized. More recently, it has been shown that 
a clear enhancement of plaquette-plaquette correlations is found by introduc- 
ing a further, third-nearest-ncighbor superexchange term J3, thus suggesting 
a possible plaquette valence-bond crystal [2"T] . 

The primary obstacle to the characterization of the phase diagram in two 
dimensions is that the lack of exact results is accompanied, in the frustrated 
case, by difficulties in applying standard stochastic numerical techniques. 
Quantum Monte Carlo methods can be applied straightforwardly only to spin- 
1/2 Hamiltonians of the form (|1.1[) . with strong restrictions on the couplings 
(e.g., J\ > and J2 < or Ji < and J2 < 0) in order to avoid a numerical 
instability known as the sign problem. This is because, in general, quantum 
Monte Carlo methods do not suffer from numerical instabilities only when it 
is possible to work with a basis in the Hilbert space where the off-diagonal 
matrix elements of the Hamiltonian are all non-positive. As an example, in a 
quantum antifcrromagnet with J\ > and J2 = on a bipartite lattice, after 
the unitary transformation 



(B being one of the two sublattices), the transformed Hamiltonian has non- 
positive off-diagonal matrix elements in the basis \x) whose states are speci- 
fied by the value of S R on each site, and J2r &r = & [22] . This implies that 
the ground state oiUTLU\ |<?o) = x ^o{ x )\ x ) , nas all-positive amplitudes, 
^(x) > 0, meaning that there exists a purely bosonic representation of the 
ground state. This property leads to the well-known Marshall-Peierls sign 
rule [221123] for the phases of the ground state of H, sign{$ a (x)} = (-l) N i< x \ 
where N^(x) is the number of up spins on one of the two sublattices. The 
Marshall-Peierls sign rule holds for the unfrustrated Heiscnbcrg model and 
even for the J1 — J2 chain at the Majumdar-Ghosh point. However, in the 
regime of strong frustration, the Marshall-Peierls sign rule is violated dramat- 
ically [24], and, because no analogous sign rule appears to exist, the ground- 
state wave function has non-trivial phases. This property turns out to be a 
crucial ingredient of frustration. 




(1.2) 
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In this respect, a very useful way to investigate the highly frustrated regime 
is to consider variational wave functions, whose accuracy can be assessed by 
employing stable (but approximate) Monte Carlo techniques such as the fixed- 
node approach |25j . Variational wave functions can be very flexible, allowing 
the description of magnetically ordered, dimerized, and spin-liquid states. In 
particular, it is possible to construct variational states with non-trivial signs 
for the investigation of the strongly frustrated regime. 

In the following, we will describe in detail the case in which the variational 
wave function is constructed by projecting fermionic mean- field states |26j . 
Variational calculations can be treated easily by using standard Monte Carlo 
techniques. This is in contrast to variational states based on a bosonic rep- 
resentation, which are very difficult to handle whenever the ground state has 
non-trivial phases [27] . Indeed, variational Monte Carlo calculations based on 
bosonic wave functions suffer from the sign problem in the presence of frustra- 
tion [28] . and stable numerical simulations can be performed only in special 
cases, for example in bipartite lattices when the valence bonds only connect 
opposite sublattices [27] . Another advantage of the fermionic representation is 
that the mean-field Hamiltonian allows one to have a simple and straightfor- 
ward representation also for the low-lying excited states (see the discussion in 
section H. 6. 1[ and also Ref. [35] for a frustrated model on a three- leg ladder). 

1.2 Symmetries of the wave function: general properties 

We define the class of projected-BCS (pBCS) wave functions on an iV-site 
lattice, starting from the ground state of a suitable translationally invariant 
BCS Hamiltonian 

Ti-BCS = _ fi8R-R') c R,a C R',<r ~ X] Ar ~ r ' C R,1 C R' ,1 + H - c - 

R,R'a R.R' 

= 5Z( £fc ~ ^ c L c fc, CT - A k c l,i c -kd + H - c -> ( L3 ) 

ka k 

where c s Ra (cr ]0 .) creates (destroys) an electron at site R with spin a, the 
bare electron band is a real and even function of k, and A k is also taken 
to be even to describe singlet electron pairing. In order to obtain a class of 
non-magnetic, translationally invariant, and singlet wave functions for spin- 
1/2 models, the ground state \BCS) of Hamiltonian l]1.3p is projected onto 
the physical Hilbcrt space of singly-occupied sites by the Gutzwiller operator 
Pg = n_R( n ^,T — n R,l) 2 ' n R,<? being the local density. Thus 

\ P BCS) = P G \BCS) = P G U(u k + v k cl A cl kd )\0), (1.4) 

k 

where the product is over all the N wave vectors in the Brillouin zone. The 
diagonalization of Hamiltonian (| 1 .3|) gives explicitly 
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I E k + e k A k I E k - ejfe 

while the BCS pairing function f k is given by 

fk = — = —ft- (1.5) 

The first feature we wish to discuss is the redundancy implied by the elec- 
tronic representation of a spin state, by which is meant the extra symmetries 
which appear when we write a spin state as the Gutzwiller projection of a 
fcrmionic state. This property is reflected in turn in the presence of a local 
gauge symmetry of the fcrmionic problem 30,31,32.. Indeed, the original spin 
Hamiltonian (|1.1[) is invariant under the local SU(2) gauge transformations 

A third transformation can be expressed in terms of the previous ones, 

S y . _^ e Wa y ^c]^ =e -i-Kcr ls /4 e i^ Xe iwa z /4 (^^j , (1.8) 

where a x , o~ y , and cr z are the Pauli matrices. As a consequence, all the different 
fermionic states connected by a local SU(2) transformation generated by (| 1 .6|) 
and (|1.7|) with site-dependent parameters give rise to the same spin state after 
Gutzwiller projection, 

P G J] EUEtn \ BCS ) = e< * P G \ BCS )' (1-9) 



where ^ is an overall phase. Clearly, the local gauge transformations defined 
previously change the BCS Hamiltonian, breaking in general the translational 
invariancc. In the following, we will restrict our considerations to the class 
of transformations which preserve the translational symmetry of the lattice 
in the BCS Hamiltonian, i.e., the subgroup of global symmetries correspond- 
ing to site-independent angles (<p,0). By applying the transformations (|1.6[) 
and (|1.7|) . the BCS Hamiltonian retains its form with modified couplings 

tR-R> t R - R > (1.10) 

A R -n, - A R . R .e 2i * (1.11) 
for Si, while the transformation £% gives 
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tR-w -> cos 26 t R _ w + i sin 6 cos 6 {A R _ R , - A* R _ RI ) 

= coa26t R ^ R > - sin26»Im/i fl _H/ (1.12) 
Ar-w -» (cos 2 QA R - R > + sin 2 6 A* R _ R ,) + i sin26t R - R > 

= RcA R _ R , +i(cos26ImA R _ R > + sm26t R _ R /) . (1.13) 

These relations are linear in t R - R / and A R - R i, and therefore hold equally 
for the Fourier components ek and A^. We note that, because A r is an even 
function, the real (imaginary) part of its Fourier transform A/~ is equal to 
the Fourier transform of the real (imaginary) part of A r . It is easy to see 
that these two transformations generate the full rotation group on the vector 
whose components are (e^, RcZifc, ImAk). As a consequence, the length of 
this vector is conserved by the full group. 

In summary, there is an infinite number of different translationally in- 
variant BCS Hamiltonians that, after projection, give the same spin state. 
Choosing a specific representation does not affect the physics of the state, 
but changes the pairing function fj~ of Eq. I|1.5[) before projection. Within 
this class of states, the only scalar under rotations is the BCS energy spec- 
trum Ek- Clearly, the projection operator will modify the excitation spectrum 
associated with the BCS wave function. Nevertheless, its invariance with re- 
spect to SU(2) transformations suggests that Ek may reflect the nature of the 
physical excitation spectrum. 

Remarkably, in one dimension it is easy to prove that such a class of wave 
functions is able to represent faithfully both the physics of Luttingcr liquids, 
appropriate for the nearest-neighbor Heiscnberg model, and the gapped spin- 
Peierls state, which is stabilized for sufficiently strong frustration. In fact, it is 
known |33| that the simple choice of nearest-neighbor hopping (e& = — 2t cos k, 
/i = 0) and vanishing gap function A/- reproduces the exact solution of the 
Haldane-Shastry model (with a gapless E^), while choosing a next-nearest 
neighbor hopping (e& = — 2t cos2fc, /i = 0) and a sizable nearest-neighbor 
pairing (Ak = A^/2t cos k) recovers the Majumdar- Ghosh state (with a gapped 

E k ) rani- 

1.3 Symmetries in the two-dimensional case 

We now specialize to the two-dimensional square lattice and investigate 
whether it is possible to exploit further the redundancy in the fcrmionic rep- 
resentation of a spin state in order to define a pairing function which breaks 
some spatial symmetry of the lattice but which, after projection, still gives a 
wave function with all of the correct quantum numbers. We will show that, 
if suitable conditions are satisfied, a fully symmetric projected BCS state is 
obtained from a BCS Hamiltonian with fewer symmetries than the original 
spin problem. For this purpose, it is convenient to introduce a set of unitary 
operators related to the symmetries of the model. 
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• Spatial symmetries: for example, lZ x (x,y) = (x,—y) and H xy (x,y) = 
(y, x). We define the transformation law of creation operators as 7Z-c x a .'R-^ 1 = 
c iz(x) a' an< ^ the action of the symmmetry operator on the vacuum is 
1Z X \0) = lZ xy \0) = |0). Note that these operators map each sublattice onto 
itself. 

• Particle-hole symmetry: PhC x a P^ 1 = i (— l) x cx,~a, where the action of 
the Ph operator on the vacuum is Ph\Q) = Yix c x ] c x JO)- 

• Gauge transformation: G c x a G^ = i c x with G |0) = |0). 

Clearly, 1Z X and lZ xy are symmetries of the physical problem (e.g., the 
Heiscnbcrg model). G is a symmetry because the physical Hamiltonian has 
a definite number of electrons, while Ph leaves invariant every configura- 
tion where each site is singly occupied if the total magnetization vanishes 
(Ni = .ATf = n). With the definition adopted, Ph acts only to multiply every 
spin state by the phase factor (—l) N i . Thus all of the operators defined above 
commute both with the Hcisenberg Hamiltonian and. because reflections do 
not interchange the two sublattices, with each other. The ground state of the 
Heisenberg model on a finite lattice, if it is unique, must be a simultaneous 
eigenstate of all the symmetry operators. We will establish the sufficient con- 
ditions which guarantee that the projected BCS state is indeed an eigenstate 
of all of these symmetries. 

Let us consider a hopping term which only connects sites in opposite sub- 
lattices, whence £k+Q = ~£fc, and a gap function with contributions from 
different symmetries (s, d x 2_ y 2, and d xy ), A = A s + A x2 ~ y2 + A xy . Further, 
we consider a case in which A s and A x ~~ v ~ couple opposite sublattices, while 
A xy is restricted to the same sublattice. In this case, the BCS Hamiltonian 

2 2 

7~t-BCS = 7~L(t, A s , A x ~ v , A xy ) transforms under the different unitary opera- 
tors acording to 

n x n(t, A s ,A x2 - y2 , a^-r- 1 = n{t, A\A x2 - y2 , -A xy ), 

K xy H(t, A% A x2 - y2 , A xy )K x l = H(t, A s ,-A x2 ~ y2 , A xy ), 
P h H(t, A s ,A x2 - y2 \ A XV )P^ 1 = H(t, A s *,A x2 - y2 \ -A xy *), 
GH{t, A s ,A x2 ~ y2 ,A xy )G- 1 = H{t, -A s , -A x2 ~ y2 7 -A xy ). 

From these transformations it is straightforward to define suitable composite 
symmetry operators which leave the BCS Hamiltonian invariant. For illustra- 
tion, in the case where A is real, one may select lZ x Ph and lZ xy if A x ~ y = 

2 2 

or lZ x Ph and 7Z xy PhG if A s = 0. It is not possible to set both A x ~ y and 
A s simultaneously different from zero and still obtain a state with all the 
symmetries of the original problem. The eigenstates \BCS) of Eq. (|1.3j) will 
in general be simultaneous eigenstates of these two composite symmetry op- 
erators with given quantum numbers, for example a x and a xy . The effect of 
projection over these states is 
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a x P G \BCS) = P G K x P h \BCS) = K x P h P G \BCS) 

= (-irilzPGlBCS), (1.14) 

where we have used that both 1Z X and P\ x commute with the projector. Anal- 
ogously, when a term A x 2_ y 2 is present, 

a xy P G \BCS) = P G Tl xy P h G\BCS) = TZ xy P h GP G \BCS) 

= K xy P G \BCS). (1.15) 

These equations show that the projected BCS state with both s and d xy , 
or d x 2_ y 2 and d xy , contributions to the gap has definite symmetry under 
reflections, in addition to being translationally invariant. The corresponding 
eigenvalues, for n = N/2 even, coincide with the eigenvalues of the modified 
symmetry operators 1Z x Ph and TZ xy PhG on the pure BCS state. 

In the previous discussion of quantum numbers, it was assumed that Uk 
and Vk arc well defined for every wave vector k. However, this condition is 
in general violated: singular fc-points are present whenever both the band 
structure e k and the gap function A k vanish, as for example with /i = 0, 
nearest-neighbor hopping and d x 2_ y 2 pairing at k = (±^,±-|). However, on 
finite lattices, this occurrence can be avoided by the choice of suitable bound- 
ary conditions. In fact we are free to impose either periodic or antiperiodic 
boundary conditions on the fermionic BCS Hamiltonian (|1.3|) . while maintain- 
ing all the symmetries of the original lattice. In our studies we have selected 
lattices and boundary conditions which do not result in singular fc-points. We 
note that the quantum numbers of the projected state do depend in general 
on the choice of boundary conditions in the fermionic BCS Hamiltonian. 

1.3.1 The Marshall-Peierls sign rule 

Another interesting property of the class of pBCS wave functions is related 
to the possibility of satisfying the Marshall-Peierls sign rule by means of a 
suitable choice of the gap function. In particular, we will restrict our con- 
siderations to the class of projected wave functions specified in Eq. (|1.4p 
when both i_R__R/ and Ar-ri are real and couple sites in opposing sublat- 
tices. We begin with the BCS Hamiltonian (| 1 .3|) and perform a particle-hole 
transformation on the down spins alone, d R ^ = c R ^ and d R ^ — e'^'^CRj, 
with Q = (jr, 7r), followed by the canonical transformation (spin rotation) 
a+(k) = (dk,-\ +*rffc,i)/v / 2 and a_(fc) = —i(dk , T -id k> i)/V2. The BCS Hamil- 
tonian then acquires the form 

H B cs = ^2[h+(k) + h-(k)}, (1.16) 

k 

where h±(k) = eka±(k)a±(k) ± iAka [ ± (k)a±(k + Q), and we have used the 
symmetry A k = —A^+q. Due to the anticommutation rules of the operators 
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a±(k), the ground state of H-bcs can be written as a tensor product of free 
states of ± fermions. Moreover, if Y,r &{ R i ' ' ' Rn)a\{Ri) ■ ■ ■ a\{R n )\Q) is the 

ground state of J2k h +{k), then J2r &*( r i ■ ■ ■ Rn)aL(Ri) ■ ■ ■ al(R„)\0) is the 
ground state of h—(k). Here we have chosen an arbitrary ordering of the 
lattice sites. The ground state of TLbcs is therefore 



J2 «W • --Rn^iK ■ ■ ■ Rn)al(Ri) ■ ■■al(R n )al(R[) ■ ■ • aL«)|0). 



R,R' 



(1.17) 

If this state is expressed in terms of the original electron operators we obtain, 
up to a factor of proportionality, 



R,R< 



&(R 1 ---R n )&*(R[---R' n ) 



ie iQRl c Rul 



iQR[ 



"1,1 



(1.18) 



In projecting over the state of fixed particle number equal to the number of 
sites, we must take the same number of creation and annihilation operators 
in the N factors of the product. The suppression of doubly occupied sites 
mandated by the Gutzwiller projector is effected by creating an up spin on 
sites where a down spin has already been annihilated. The only terms which 
survive are then those with {R} = {R 1 }, namely 



R 



\9(Rl---Rn)\ 



e c R 1 ,1 c Ri,l 



"i?„,T Ci W ' ' ' c i,J. 



J 



|0). 



(1.19) 

Finally, by moving the down-spin creation operators to the left, one may order 
the operators according to the specified ordering of the sites in the lattice, 
independently of the spin, without introducing any further phase factors. On 
this basis, the resulting wave function has exactly the Marshall-Peierls sign. 



1.3.2 Spin correlations 



Finally, we would like to calculate the form of the long-range decay of the 
spin correlations in a BCS state. Here, we will show only that the pure BCS 
state before projection is characterized by correlations which maintain the 
symmetries of the lattice even when the BCS Hamiltonian breaks the reflection 
symmetries due to the presence of both A x ~ v and A xy couplings. Because 
the BCS state (|1.4[) is a translationally invariant singlet, it is sufficient to 
calculate the longitudinal correlations (S^Sq). A straightforward application 
of Wick's theorem leads (for r ^ 0) to 

(S z r S%) K-[g 2 (r) + h 2 (r)] , (1.20) 

g(r) = Jd 2 k^-e lkr , (1.21) 
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h(r) = f d 2 k — e lkr . (1.22) 
J E k 

Note that when the gap function A k has both d x 2_ y 2 and d xy contribu- 
tions, the correlation function apparently breaks rotational invariance. Equa- 
tion (|1.20p can be written equivalently in Fourier space as 

S(q) = (S*S*_ a ) oc I d 2 k , e ^ k+q + A k A k+q (J ^ 

f 2 + A 2 

Now the effect of an cc-reflection 1Z X on the wave vector q can be deduced by 
setting A k = A\ ~ v +A 3 k yy and changing the dummy integration variable k — > 

Tl x (k + Q), whence A k+Q = -Af'^+A^ and A R „ k = Z\f ^ '- Af '. The 
net result of these transformations is simply S(7Z x q) = S(q), demonstrating 
that the spin correlations of a BCS state are isotropic, even if the gap function 
breaks rotational invariance before Gutzwiller projection. 

The explicit evaluation of the long-range decay of g{r) for a d x 2_ y 2 gap 
shows that spin correlations in a BCS state (i.e., before projection) display a 
power-law decay due to the presence of gapless modes: (S z Sg) ~ 1/r 4 for sites 
on opposite sublattices, while (S*S§) vanishes for sites on the same sublattice. 
A similar result is also expected in the presence of a finite A^ v , because gapless 
modes are present also in this case. 



1.4 Connection with the bosonic representation 

We turn now to a detailed discussion of the relation between the fermionic [55] 
and bosonic [57] representations of the RVB wave function. Recently, bosonic 
RVB wave functions have been reconsidered by Beach and Sandvik [Ml [551 
[36j[37]. In particular, it has been possible to improve the earlier results of 
Ref. |27j . either by assuming some asymptotic form of the bond distribu- 
tion |37j or by unconstrained numerical methods |34| . This wave function has 
been demonstrated to be extremely accurate for the unfrustrated model with 

j 2 = o [HIM]. 

In the fermionic representation, we have 



\pBCS) = P G exp 



/B.*'(4,T c fl', 



R<R' 



|0>, (1.24) 



where Pq projects onto the physical subspace with one electron per site and 
fR,R> is the pairing function, given by the Fourier transform of Eq. (| 1 . 5|) . 
The constraint R < R' implies the definition of an (arbitrary) ordering of 
the lattice sites: here and in the following, we will refer to the lexicographical 
order. For simplicity, let us denote the singlet operator as &n.,n' = (cr +c' r , ^ + 
C R' \ C R i) - Once the Gutzwiller projector is taken into account, we have that 
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\ P BCS) = Y, E fR i< R 'i ■ ■ ■ fRn,K e RuR[ ■ ■ ■ ®Rn,R' n |0>, (1-25) 

i?i<...<i?„ P(R') 

where n — N/2 and P(R') represents the permutations of the n sites R' k not 
belonging to the set {R}, satisfying Rk < R' k for every k. The sum defines all 
the (N — 1)!! partitions of the N sites into pairs. 

On the other hand, the bosonic RVB wave function may be expressed in 
terms of the spin-lowering operator, S R , as 

\RVB)= Yl E fRuR[---fR:,K(^ 1 -SnJ...(S^-S^ n )\F), 

R 1 <...<R n P(R') 

(1.26) 

where the sum has the same restrictions as before and \F) is the (fully po- 
larized) ferromagnetic state. In the bosonic representation, a valence-bond 
singlet is antisymmetric on interchanging the two sites and, therefore, a direc- 
tion must be specified. The condition Rk < R' k fixes the phase (i.e., the sign) 
of the RVB wave function. In order to make contact between the two repre- 
sentations, we express \F) and S R in terms of fermionic operators, namely 
\F) = c| >t . . -cjy^lO) and S R = c Rl c R .p Then 

Ri<...<R n P(R') 

(1.27) 

where £{r.r'} = ±1 is a configuration-dependent sign arising from the re- 
ordering of the fermionic operators (1, . . . N) —* (R\, R[ . . . R n , R' n ). The two 
representations are therefore equivalent only if 

£{R,R'} /fli.fii ■ • • fR n ,R!n = fRT,R[ ■ ■ ■ flilR' n (!- 28 ) 

for all the valence-bond configurations. In general, for a given f R os R , , this con- 
dition cannot be satisfied by any choice of Jr^r/. Remarkably, this is however 
possible for the short-range RVB state ;38,.39j, where only nearest- neighbor 
sites are coupled by f R °^ R > = 1- Indeed, by using the Kasteleyn theorems [40] . 
it is possible to prove that Eq. (|1.28p can be fulfilled on all planar graphs (for 
example in short-range RVB states on lattices with open boundary condi- 
tions). In fact, the left-hand side of Eq. (|1.28p is a generic term in the Pfaffian 
of the matrix 

M(R,R') = S f f M ' ^^2' (!- 29 ) 
\ -fw,R for R> R- { ' 

As a consequence, following the arguments of Kasteleyn, it is always possible 
to orient all the bonds in such a way that in all cycles of the transition graph 
the number of bonds oriented in either directions is odd |40j . Notice that the 
latter way to orient the bonds will in general be different from the one used 
in Eq. (|1.26p . Thus we define = 1 (with R < R') if the bond is oriented 

from R to R' , and {r^ri = — 1 otherwise. In summary, in order to define the 
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fcrmionic pairing function fnw once we know the oriented planar graph, it 
is necessary to: 

• label the sites according to their lexicographical order, 

• orient the bonds in order to meet the Kasteleyn prescription, and 

• take fn t R> = 1 for the bond oriented from R to i?', and fB„R/ = — 1 
otherwise. 

This construction is strictly valid only for planar graphs, namely for graphs 
without intersecting singlets, implying that open boundary conditions must 
be taken. In this case it is known that a unique short-range RVB state can 
be constructed. Periodic boundary conditions imply the existence of four de- 
generate states, which are obtained by inserting a cut (changing the sign of 
the pairing function on all bonds intersected) that wraps once around the 
system, in the x, y or both directions [39]. These different states have the 
same bulk properties and, despite the fact that it would be possible to obtain 
a precise correspondence between bosonic and fermionic states, their physical 
properties can be obtained by considering a single (bosonic or fermionic) wave 
function. 



1.5 Antiferromagnetic order 

In the preceding sections we have considered the mean-field Hamiltonian ()1.3j) 
containing only hopping and pairing terms. In this case, even by considering 
the local SU(2) symmetries described above, it is not possible to generate a 
magnetic order parameter. The most natural way to introduce an antiferro- 
magnetic order is by adding to the BCS Hamiltonian of Eq. p. 31) a magnetic 
field 

T^bcs+af — 7~Ibcs + 7~Laf- (1.30) 

Usually, the antiferromagnetic mean-field order parameter is chosen is chosen 
to lie along the z-direction [41], 

Haf = A AF ^ e^-«(4 iT CR )T - cJ^cju), (1.31) 

R 

where Q is the antiferromagnetic wave vector (e.g., Q = for the Neel 

state). However, in this case, the Gutzwiller-projected wave function obtained 
from the ground state of Eq. (|1.30p overestimates the correct magnetic order 
parameter (see section lT.6.2p . because important quantum fluctuations are ne- 
glected. A more appropriate description which serves to mitigate this problem 
is obtained by the introduction of a spin Jastrow factor J which generates 
fluctuations in the direction orthogonal to that of the mean-field order pa- 
rameter [4"2ll4"3"] . Therefore, we take a staggered magnetic field Aaf along the 
x axis, 
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U AF = A AF e^- R (c^ t c Rd + J Rd c RA ), (1.32) 

R 

and consider a long-range spin Jastrow factor J 

^ = cx P^E^-«'^'j ' (1-33) 

vr-ri being variational parameters to be optimized by minimizing the energy. 
The Jastrow term is very simple to compute by employing a random walk 
in the configuration space \x) = c R . . . ct a |0) defined by the electron 
positions and their spin components along the z quantization axis, because 
it represents only a classical weight acting on the configuration. Finally, the 
variational ansatz is given by 

\ P BCS + AF) = JP Sz =oP G \BCS + AF), (1.34) 

where Ps^=o is the projector onto the S z = sector and \BCS + AF) is 
the ground state of the Hamiltonian (| 1 .30[) . It should be emphasized that this 
wave function breaks the spin symmetry, and thus, like a magnetically ordered 
state, it is not a singlet. Nevertheless, after projection onto the subspace with 
S% ot = 0, the wave function has (S R ) = (S R ) = (S R ) = 0. Furthermore, the 
correlation functions (S R S RI ) and (S R S RI ) have the same behavior, and hence 
the staggered magnetization lies in the x—y plane [32] . 

The mean-field Hamiltonian (|1.30p is quadratic in the fermionic operators 
and can be diagonalized readily in real space. Its ground state has the general 
form 

\BCS + AF) = cxp I \ fRiAA',.- I |0>, (1-35) 

\ R,R',a.cr' J 

where the pairing function f RR , is an antisymmetric 2N x 2N matrix. We 
note that in the case of the standard BCS Hamiltonian, with Aaf = or even 
with Aaf along z, f RR , = f R R i = 0, whereas in the presence of a magnetic 
field in the x~y plane the pairing function acquires non-zero contributions 
also in this triplet channel. The technical difficulty when dealing with such a 
state is that, given a generic configuration with a definite z-component of the 
spin, \x) = c Ruai . ■■c Rn (7n \0), one has 

(x\BCS + AF) = Pf[F], (1.36) 

where Pf[F] is the Pfaffian of the pairing function 



/(Jo R'a\ T> -ft/?) 



(1.37) 



in which the matrix F has been written in terms of N x N blocks and R a and 
R' a are respectively the positions of the up and down spins in the configuration 
\x) Hfl. 
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Fig. 1.1. Left panels: comparison between exact (empty circles) and variational (full 
dots) results for the spin structure factor S(q) on a chain with 30 sites. Right panel: 
variational results for the spin structure factor S(q) for 122 sites and J2/J1 > 0.5. 
Inset: position of the maximum of S(q), indicated by 9, as a function of the ratio 
J2/J1 (full dots). For comparison, the DMRG results of Ref. [8] are also shown 
(empty triangles). 

1.6 Numerical Results 

In this section, we report numerical results obtained by the variational Monte 
Carlo method for the one- and two-dimensional lattices. The variational pa- 
rameters contained in the BCS and BCS+AF mean-field Hamiltonians of 
Eqs. (|1.3|) and (|1.30|) . as well as the ones contained in the spin Jastrow 
factor, (|1.33p . can be obtained by the optimization technique described in 
Ref. 05]. 

1.6.1 One-dimensional lattice 

We begin by considering the one-dimensional case, where the high level of ac- 
curacy of the pBCS wave function can be verified by comparison with Lanczos 
and DMRG results. We consider the Hamiltonian (jl.ip on a chain with N 
sites and periodic boundary conditions, and first discuss in some detail the 
parametrization of the wave function. For J 2 = 0, a very good variational state 
is obtained simply by projecting the free-electron Slater determinant, where 
€k = —2tcosk |46| . Then, in one dimension, the nearest-neighbor BCS pair- 
ing A\ is irrelevant, and, in order to improve the variational energy, a third- 
neighbor BCS pairing A3 must be considered in addition; a second-neighbor 
pairing term A2, like the chemical potential, violates the Marshall-Peierls sign 
rule, which must hold at J2 = 0, and thus is not considered. To give some 
indication of the accuracy of the wave function, we note that for N = 30, 
the energy per site of the projected Fermi sea is E/Ji = —0.443060(5), while 
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by optimizing A\ and A 3 one obtains Ej J\ = —0.443934(5), the exact result 
being £ /J x = -0.444065. 

When both the chemical potential and A 2 vanish, the particle-hole trans- 
formation Ph (section ll.3[) is a symmetry of the BCS Hamiltonian. In finite 
chains, the BCS ground state is unique only if the appropriate boundary con- 
ditions are adopted in Hbcs- if, for example, N = 41 + 2 with integer 
periodic boundary conditions (PBC) should be used, while the imposition of 
antiperiodic boundary conditions (APBC) causes four zero-energy modes to 
appear in the single-particle spectrum. By filling these energy levels, we can 
form six orthogonal BCS ground states in the S z = subspace, which, in the 
thermodynamic limit, are degenerate with the ground state of the BCS Hamil- 
tonian with PBC. However, two of these states have the wrong particle-hole 
quantum number and are therefore annihilated by the Gutzwiller projector. 
If the remaining four BCS states (three singlets and one triplet) arc still 
orthogonal after projection, one may infer either the presence of a gapless 
excitation spectrum or of a ground-state degeneracy. We have built these five 
projected states (one with PBC and four with APBC) for a N ~ 30 chain 
and variational parameters appropriate for J 2 = 0. Two of them belong to the 
symmetry subspace of the ground state and represent the same physical wave 
function (their overlap is Ktfq l^)! = 0.999), two of them are singlets with mo- 
mentum it relative to ground state and again show an extremely large overlap 
(K^l^)! = 0.921), and the remaining state is a triplet with momentum tt 
relative to the ground state. Therefore only three independent states can be 
obtained by this procedure. It is remarkable that by optimizing the parame- 
ters for the ground state, and without any additional adjustable parameters, 
these three variational states have overlap higher than 98.7% with the exact 
cigenstates of the Hciscnbcrg Hamiltonian in the lowest levels of the confor- 
mal tower of states [47] , thereby reproducing with high accuracy the ground 
state and the lowest singlet and triplet modes. 

By increasing the frustrating interaction, the parameters A\ and A3 (both 
real) grow until a divergence at J2/J1 ~ 0.15. For larger values of J2/J1, the 
band structure changes: here eu = — 2i'cos2fc — [i and a non-vanishing BCS 
pairing is found, leading to a finite gap in the BCS spectrum, Ek = \J ej, + A\ . 
In this regime, although the variational wave function is translationally in- 
variant, it shows a long-range order in the dimer-dimer correlations (see be- 
low) . Similar behavior has been also discussed in Ref. [I5J for a complex wave 
function on ladders with an odd number of legs. The variational parameters 
appropriate for this regime correspond to a gapped BCS single-particle spec- 
trum for both PBC and APBC, and then only two states can be constructed. 
However, the symmetry subspace of the variational wave function depends on 
the choice of boundary conditions, implying a ground-state degeneracy. In a 
chain of N = 30 sites and for J2/J1 = 0.4, we found that the two singlets 
which collapse in the thermodynamic limit (due to the broken translational 
symmetry) have overlaps higher than 99% with the two pBCS wave functions 
corresponding to the same variational parameters and different boundary con- 
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Fig. 1.2. Left panels: dimer-dimer correlations as a function of distance for exact 
(empty circles) and variational (full dots) calculations on a chain of 30 sites. Right 
panel: dimer order parameter of Eq. (|1.40p as a function of the ratio J2/J1 for 
N — 30, 50, and 150; the extrapolation to the N — > 00 limit is also shown, together 
with the DMRG results of Ref. [8]. 



ditions. This shows that the pBCS class of wave functions is able to describe 
valence-bond crystals and broken-symmetry states. By increasing further the 
ratio J2/J1, beyond 0.5 wc found that, while the bare dispersion is again 
efc = — 2t cosfc, Ai acquires a finite value (together with A\ and ^3), showing 
both dimerization and short-range incommensurate spin correlations. 

The primary drawback of the variational scenario is that the critical point 
for the transition from the gapless fluid to the dimerized state is predicted 
around J2/J1 ~ 0.15 (where the best singlet variational state is the fully 
projected Fermi sea), considerably smaller than the known critical point 
J2/J1 ~ 0.241. This estimate does not change appreciably on considering 
further parameters in the BCS Hamiltonian (jl.3|) . probably because the vari- 
ational wave function does not describe adequately the backscattering term 
which is responsible for the transition |49j . In order to improve this aspect, 
it is necessary to include the spin Jastrow factor of Eq. (|1.33j) (without the 
mean- field magnetic parameter Aaf)- In this way, although the variational 
state is no longer a singlet, the value of the square of the total spin (S 2 ) re- 
mains very small (less than 0.002 and 0.02 for 30 and 122 sites, respectively) 
and no long-range magnetic order is generated. The Jastrow factor is particu- 
larly important in the gapless regime: despite the fact that the gain in energy 
with respect to the singlet state is less than 10 -4 J\ (specifically, for J2 = wc 
obtain EjJ\ = —0.444010(5)), this correction is able to shift the transition, 
always marked by the divergence of the BCS pairings, to J2/J1 ~ 0.21, a 
value much closer to the exact result. A finite value of the chemical potential 
is generated for 0.22 < J2/J1 < 0.5. 



1 Variational wave functions for frustrated magnetic models 17 

Let us now investigate the physical properties of the variational wave func- 
tion by evaluating some relevant correlation functions. The spin structure 
factor is defined as 

s («)4^ e<5(K " B,)< ™- (i - 38) 

R,R' 

While true long-range magnetic order cannot be established in one-dimensional 
systems, for J2/J1 <C 1 the ground state is quasi-ordered, by which is meant 
that it sustains zero-energy excitations and S(q) displays a logarithmic diver- 
gence at q = it. In Fig. we show the comparison of the spin structure fac- 
tor for an exact calculation on TV = 30 and for the variational wave function. 
Remarkably, the variational results deliver a very good description of S(q) 
in all the different regimes: for small J2/J1, where the spin fluctuations are 
commensurate and there is a quasi-long-range order, for 0.21 < J%j J\ < 0.5, 
where the spin fluctuations are still commensurate but short-range, and for 
Jil J\ > 0.5, where they are incommensurate and the maximum of S(q) moves 
from q = ir at J2/ J\ = 0.5 to q = tt/2 for J2/J1 — ► 00. Indeed, it is known that 
the quantum case is rather different from its classical counterpart [5]: while 
the latter shows a spiral state for J2/J1 > 0.25, with a pitch angle 9 given by 
cosf? = — J1/4J2, the former maintains commensurate fluctuations at least up 
to the Majumdar-Ghosh point. The behavior of S(q) for a large lattice with 
122 sites and J2/J1 > 0.5 is shown in Fig. 11.1) where we find good agreement 
with previous numerical results based upon the DMRG technique [5] . 

In the one-dimensional J1 — J2 model, there is clear evidence for a Berezinskii- 
Kosterlitz-Thouless transition on increasing the ratio J2/J1 from a gapless 
Luttingcr liquid to a dimerized state that breaks the translational symme- 
try. In order to investigate the possible occurrence of a dimerized phase, we 
analyze the dimer-dimer correlation functions of the ground state, 

0(R — R') = (S^S R+X S^,S^, +X ) — (S R Sx +x )(Sft,Sft, +x ). (1.39) 

While this definition considers only the z component of the spin operators, in 
the presence of a broken spatial symmetry the transverse components must 
also remain finite at large distances, displaying also a characteristic alter- 
nation. By contrast, in the gapless regime, the dimer correlations decay to 
zero at large distances. The differing behavior of these correlations is easy 
to recognize, with oscillatory power-law decay in the Luttingcr regime and 
constant-amplitude oscillations in the dimerized phase. Figure fl~2l illustrates 
the comparison of the dimer-dimer correlations (|1.39p between the exact and 
the variational results on a chain with 30 sites. Also for this quantity we obtain 
very good agreement for all values of the frustrating superexchange J2, both 
in the gapless and in the dimerized regions. Following Ref. [8], it is possible 
by finite-size scaling to obtain an estimate of the dimer order parameter from 
the long-distance behavior of the dimer-dimer correlations, 

d 2 = 9 lim \(O(R-x)-20(R)+G(R + x)\, (1.40) 

|-R|— »oo 
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Fig. 1.3. Upper left panel: energy per site as a function of cluster size N, showing 
exact results (full squares), variational results obtained by considering Eqs. (|1.30|l 
and (|1.32[) with the spin Jastrow factor (|1.33[) (full circles), and variational results 
obtained with Eqs. (|1.30|l and (|1.31[) (full triangles). The results obtained by optimiz- 
ing the pairing function f^ s R i of the bosonic representation, described in section [L4l 
are also shown (empty circles) [53]. Lower left panel: staggered magnetization with 
the same notation as in the upper panel. 

Right panel: static structure factor S(q) for a cluster with N — 242 (tilted by 45°): 
variational results for the state of Eqs. (|l,30[l and (|1.32[) with a long-range Jastrow 
factor (full dots) and for the wave function of Eqs. (|1,30|1 and (|1.31[) with a nearest- 
neighbor Jastrow factor (empty triangles). Lower inset: detail at small momenta. 
Upper inset: square of total spin (S 2 ) as a function of N for the two states, using 
the same symbols. 



where the factor 9 is required to take into account the fact that in Eq. (jl.39[) wc 
considered only the z component of the spin operators. In Fig. 11.21 we present 
the values of the dimer order parameter as a function of J2/J1 for three 
different sizes of the chain, and also the extrapolation in the thermodynamic 
limit, where the agreement with the DMRG results of Rcf. [5] is remarkable. 

1.6.2 Two-dimensional lattice 

We move now to consider the two-dimensional case, starting with the unfrus- 
trated model (J2 = 0), for which exact results can be obtained by Monte 
Carlo methods [SDl[nU[52]. In the thermodynamic limit, the ground state 
is antiferromagnetically ordered with a staggered magnetization reduced to 
approximately 60% of its classical value, namely M ~ 0.307 j?Tl[5^] . This 
quantity can be obtained both from the spin-spin correlations at the largest 
distances and from the spin structure factor S(q) at q = (n,ir). In the fol- 
lowing, we will consider the former definition and will calculate the isotropic 
correlations (S^ • S^'), because this quantity is known to have smaller finite- 
size effects [SHIED]- For the unfrustrated case, the best wave function has 
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Fig. 1.4. Left panel: spin-spin correlations at the largest distances as a function of 
the ratio J2/J1 for different cluster sizes N. Inset: finite-size scaling for J2/J1 = 0.42 
(empty triangles), 0.44 (full triangles), 0.46 (empty squares), 0.48 (full squares), 
0.50 (empty circles), and 0.52 (full circles). Right panel: spin-spin correlations at 
the largest distances for J2/J1 > 0.7. 



gfc = — 2t(cosk x + cosky) and a pairing function with d x 2_ y 2 symmetry, 
2 2 

~ V = ^1 ( cos ~ cos ky) (possibly also with higher harmonics connecting 
opposite sublatticcs) . The quantity Aj^p in Eq. (|1.32p has a finite value and 
the spin Jastrow factor (|1.33p has an important role. 

Figure 11.31 shows the comparison of the variational calculations with the 
exact results, which are available for rather large system sizes. In the un- 
frustrated case, the bosonic representation is considerably better than the 
fermionic one: the accuracy in the energy is around 0.06% and the sublat- 
tice magnetization is also very close to the exact value [341E3]. However, the 
fermionic state defined by Eqs. (| 1 .30)) and p.32[) . in combination with the spin 
Jastrow factor, also provides a very good approximation to the exact results 
(energy per site and staggered magnetization), whereas the wave function de- 
fined by Eqs. Q1.30P and (|1.31|1 is rather inaccurate. It should be emphasized 
that when the Jastrow factor is included, the slopes of the finite-size scaling 
functions are also remarkably similar to the exact ones, both for the energy 
per site eo and for the magnetization M. This implies that the pBCS wave 
function provides an accurate estimate of the spin velocity c, of the transverse 
susceptibility x_i_, and as a consequence of the spin stiffness, p s = c 2 x±- By 
contrast, the wave function without the Jastrow factor leads to a vanishing 
spin velocity. We note that in this case the staggered magnetization M ~ 0.365 
is also overestimated in the thermodynamic limit. 

The functional form of the Jastrow factor at long ranges, which can be 
obtained by minimizing the energy, is necessary to reproduce correctly the 
small-g behavior of the spin-structure factor S(q), mimicking the Goldstone 
modes typical of a broken continuous symmetry |42| . Indeed, it is clear from 
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Fig. ll.3l that only with a long-range spin Jastrow factor it is possible to obtain 
S(q) ~ |g| for small momenta, consistent with a gapless spin spectrum. By 
contrast, with a short-range spin Jastrow factor (for example with a nearest- 
neighbor term), S(q) ~ const for small q, which is clearly not correct [12]. 
Finally, it should be emphasized that the combined effects of the magnetic 
order parameter Aaf and the spin Jastrow factor give rise to an almost 
singlet wave function, strongly reducing the value of (S 2 ) compared to the 
case without a long-range Jastrow term (see Fig. II. 3j) . 

On increasing the value of the frustrating superexchange J2, the Monte 
Carlo method is no longer numerically exact because of the sign problem, 
whereas the variational approach remains easy to apply. In Fig. II. 41 we present 
the results for the spin-spin correlations at the maximum accessible distances 
for J2/J1 < 0.52. It is interesting to note that when J2/J1 > 0.4, a sizable 
energy gain may be obtained by adding a finite pairing connecting pairs on 

the same sublattice with d xy symmetry, namely Ak = A\ ~ v + A^ y [54]. The 
mean-field order parameter Aaf remains finite up to J2/J1 ~ 0.5, whereas 
for J2/J1 > 0.5 it goes to zero in the thermodynamic limit. Because the 
Jastrow factor is not expected to destroy the long-range magnetic order, the 
variational technique predicts that antiferromagnetism survives up to higher 
frustration ratios than expected [12] , similar to the outcome of a Schwinger bo- 
son calculation |55j . The magnetization also remains finite, albeit very small, 
up to J2/J1 = 0.5 (see Fig. I1.4[) . We remark here that by using the bosonic 
RVB state, Beach argued that the Marshall-Peierls sign rule may hold over a 
rather large range of frustration, namely up to J2/J1 = 0.418, also implying 
a finite staggered magnetic moment |37| . In this approach, if one assumes a 
continuous transition from the ordered to the disordered phase, the critical 
value is found to be J2/J1 = 0.447, larger than the value of Ref. [T2] and much 
closer to our variational prediction. We note in this context that recent results 
obtained by coupled cluster methods are also similar, i.e., J2/J1 ~ 0.45 for 
a continuous phase transition between a Neel ordered state and a quantum 
paramagnet [56j . 

In the regime of large J2/J1 (i.e., J2/J1 > 0.65), collinear order with pitch 
vectors Q = (0, tt) and Q = (tt, 0) is expected. The pBCS wave function 
is also able to describe this phase through a different choice for the bare 
electron dispersion, namely = — 2t'[cos(k x + k y ) + cos(fc x — k y )] and Ak = 
Aicosk x + Zi2[cos(fc x + k y ) — cos(k x — k y )], with A\ — * for J1/J2 — ► 0. 
Further, the antifcrromagnetic wave vector Q in Eq. (|1.32p is Q = (n, 0). The 
variational wave function breaks the reflection symmetry of the lattice and, 
in finite systems, its energy can be lowered by projecting the state onto a 
subspacc of definite symmetry. The results for the spin-spin correlations are 
shown in Fig. 11.41 By decreasing the value of J2/J1, we find clear evidence of 
a first-order phase transition, in agreement with previous calculations using 
different approaches [TTIUS] ■ 

For 0.5 < J2/J1 < 0.65, the best variational wave function has no mag- 
netic order (Aaf = and no Jastrow factor) and the BCS Hamiltonian has 
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Fig. 1.5. Comparison between exact and variational results on a 6 x 6 lattice. 
The pBCS wave function has Aaf = and no Jastrow factor. Upper panel: 
average sign (s) of Eq. (|l,41[l (full circles); the Marshall- Peierls sign, (s)ma = 
J2 X I (a: I *b) | 2 sign { (a;|<f r o)(-l) JVT(a:> } , is also shown (full triangles). Middle panel: ac- 
curacy of the ground-state energy, AE/Eq = {Eq — E p bcs) I 'Eq, where Eq and 
EpBCS are the exact and the variational energies, respectively. Lower panel: overlap 
between the exact \\Po) and variational \pBCS) states (full circles). The norm of 
the projection of the variational state onto the subspace spanned by the two lowest- 
energy states in the same symmetry sector is also shown (full squares) close to the 
first-order transition to the collinear state. 

gfc = — 2t{cosk x + cosky) and A}. = A x k ~ v + A x k v , where A k ~ v connects 
pairs on opposite sublattices while A x k v is for same sublatticc. With this spe- 
cific electron pairing, the signs of the wave function are different from those 



Table 1.1. Energies per site for a 6 x 6 lattice. E v bcs obtained from the pBCS 
wave function (with Aaf = and no Jastrow factor), Elr-rvb from the long- 
range bosonic RVB state, optimizing just one parameter using the master-equation 
method [58], and Esr-rvb obtained by diagonalizing the J1 — J2 model in the 
nearest-neighbor valence-bond basis [57]. The exact results Eq are also reported. 



J2/J1 


Esr-rvb 


Elr-rvb 


EpBCS 


Eq 


0.30 


-0.54982 


-0.5629(5) 


-0.55569(2) 


-0.56246 


0.35 


-0.53615 


-0.5454(5) 


-0.54134(1) 


-0.54548 


0.40 


-0.52261 


-0.5289(5) 


-0.52717(1) 


-0.52974 


0.45 


-0.50927 




-0.51365(1) 


-0.51566 


0.50 


-0.49622 




-0.50107(1) 


-0.50381 


0.55 


-0.48364 




-0.48991(1) 


-0.49518 


0.60 


-0.47191 




-0.47983(2) 


-0.49324 
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Fig. 1.6. Left panels: dimer-dimer correlations as a function of the Manhattan 
distance for exact (empty circles) and variational (full dots) calculations on a 6 x 6 
cluster. Right panel: finite-size scaling of the dimer order parameter for J2/J1 = 0.55. 



predicted by the Marshall-Peierls rule and are much more similar to the exact 
ones. We define 

(s) =Y / \(x\pBCS)\ 2 sign{(x\ P BCS)(x\%)}, (1.41) 

X 

where \pBCS) and are the variational and the exact states, respectively. 
This quantity is shown in Fig. 11.51 together with the Marshall-Peierls sign, for 
a 6 x 6 lattice. The variational energy, the very large overlap with the exact 
ground state, and the dimer-dimer correlations shown in Figs. 11.51 and 11.61 
all reflect the extremely high accuracy of this state in the strongly frustrated 
regime. On small clusters, the overlap between the variational wave function 
and the ground state deteriorates for J2/J1 > 0.55. This may be a conse- 
quence of the proximity to the first order transition, which marks the onset 
of collinear magnetic order, and implies a mixing of the two finite-size ground 
states corresponding to the coexisting phases. 

In Table II. 1( we report the comparison between the energies of the non- 
magnetic pBCS wave function and two bosonic RVB states. The first is ob- 
tained by a full diagonalization of the J1 — J2 model in the nearest-neighbor 
valence-bond basis, namely by optimizing all the amplitudes of the indepen- 
dent valence-bond configurations without assuming the particular factorized 
form of Eq. (|1.26j) [57] , Although this wave function contains a very large 
number of free parameters, its energy is always higher than that obtained 
from the pBCS state, showing the importance of having long-range valence 
bonds. A further drawback of this approach is that it is not possible to perform 
calculations on large system sizes, the upper limit being N ~ 40. The second 
RVB state is obtained by considering long-range valence bonds, with their 
amplitudes given by Eq. (|1 .26|) and optimized by using the master-equation 
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Fig. 1.7. Upper panel: phase diagram of the J1—J2 model on the square lat- 
tice, as deduced from the variational approach. Lower panel: phase diagram of the 
anisotropic triangular lattice from Ref. [64]. The approximate locations of some 
relevant materials are indicated by the arrows. 



scheme [58] . While this wave function is almost exact in the weakly frustrated 
regime, its accuracy deteriorates on raising the frustrating interaction, and 
for J2/J1 > 0.425 the minus-sign problem precludes the possibility of reliable 
results. On the other hand, the pBCS state (without antiferromagnetic order 
or the Jastrow term) becomes more and more accurate on approaching the 
disordered region. Remarkably, for J2/J1 = 0.4, the energy per site in the 
thermodynamic limit obtained with the long-range bosonic wave function is 
E/Ji = —0.5208(2), which is very close to and only slightly higher than that 
obtained from the fcrmionic representation, E/Ji = —0.5219(1). 

In the disordered phase, the pBCS wave function does not break any lattice 
symmetries (section [L3]| and does not show any tendency towards a dimer- 
ization. Indeed, the dimer order parameter d (calculated from the correlations 
at the longest distances) vanishes in the thermodynamic limit, as shown in 
Fig. 11.61 implying a true spin-liquid phase in this regime of frustration. This 
fact is in agreement with DMRG calculations on ladders with odd numbers 
of legs, suggesting a vanishing spin gap for all values of J2/J1 [59], in sharp 
contrast to the dimerized phase, which has a finite triplet gap. 

Taking together all of the above results, it is possible to draw the (zero- 
temperature) phase diagram generated by the variational approach, and this 
is shown in Fig. 11.71 

We conclude by considering the important issue of the low-energy spec- 
trum. In two dimensions, it has been argued that the ground state of a spin- 1/2 
system is cither degenerate or it sustains gapless excitations [60], in analogy 
to the one-dimensional case [61] . In Ref. [62] , it has been shown that the wave 
function with both d x 2 _ y 2 and d xy parameters could have topological order. In 
fact, by changing the boundary conditions of the BCS Hamiltonian, it should 
be possible to obtain four different projected states which in the thermody- 
namic limit are degenerate and orthogonal but, however, not connected by 
any local spin operator. In this respect, it has been argued more recently that 



24 



Federico Becca, Luca Capriotti, Alberto Parola, and Sandro Sorella 




Fig. 1.8. Nearest-neighbor pairing function consistent with the sign convention of 
the short-range RVB state in the triangular lattice: solid (dashed) lines represent 
positive (negative) values. Note that the unit cell contains two sites, indicated by 
empty and full circles. 

a topological degeneracy may be related to the signs of the wave function and 
cannot be obtained for states satisfying the Marshall-Peierls rule [S3] ■ 

In the spin-liquid regime, the simultaneous presence of A^ ~ y and A k 
could shift the gaplcss modes of the unprojccted BCS spectrum from 
(±7r/2, ±7r/2) to incommensurate fc-points along the Fermi surface determined 
by efe = 0. However, we have demonstrated recently that a particular Af v pair- 
ing, A^ v oc sin(2fc x ) sin(2fc y ), may be imposed, in order to fix the nodes at 
the commensurate points (±7r/2, ±7t/2), without paying an additional energy 
penalty. Once Ek is connected to the true spin excitations, a gapless spec- 
trum is also expected. At present, within a pure variational technique, it is 
not possible to assess the possibility of incommensurate, gapless spin excita- 
tions being present. An even more challenging problem is to understand if the 
topological states could survive at all in the presence of a gapless spectrum. 

1.7 Other frustrated lattices 

In this last section, we provide a brief overview of related variational studies 
performed for other lattice structures. In particular, wc discuss in some detail 
the symmetries of the variational wave function on the anisotropic triangu- 
lar lattice, considered in Ref. [64] . In this case, one-dimensional chains with 
antiferromagnetic interaction J are coupled together by a superexchange J', 
such that by varying the ratio J' /J, the system interpolates between decou- 
pled chains (J 7 = 0) and the isotropic triangular lattice (J 7 = J); the square 
lattice can also be described in the limit of J = 0. The case with J' < J may 
be relevant for describing the low-temperature behavior of CS2CUCI4 [65], 
whereas J' ~ J may be pertinent to the insulating regime of some organic 
materials, such as n — (ET)2Cu2(CN)3 [66] . 

In Ref. [64], it has been shown that very accurate variational wave func- 
tions can be constructed, providing evidence in favor of two different spin- 
liquid phases, a gapped one close to the isotropic point and a gapless one 
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Fig. 1.9. Left panel: signs of the real hopping terms Xr.r' of the U(l) Dirac spin 
liquid on the kagome lattice [Eq. (|l,43|l ] 69 . Solid (dashed) lines represent positive 
(negative) values. The unit cell contains six sites (three empty and three full circles 
inside the boxes) . Right panel: nearest- neighbor pairing function consistent with the 
sign convention of the short-range RVB state in the kagome lattice: solid (dashed) 
lines represent positive (negative) values. The unit cell also contains six sites in this 
case. 

close to the one-dimensional regime, see Fig. 11.71 We focus our attention on 
the isotropic point. In this case, a natural variational ansatz is the bosonic 
short-range RVB state of Eq. (|1.26[) [2~5] . Exact numerical calculations for the 
6x6 isotropic model have shown that the overlap between the short-range 
RVB wave function and the ground state is very large, | (RVB\^o) | 2 = 0.891, 
and also that the average sign (s) = J2x\( x \'^ r o)\ 2s ' 1 S n {( x \^ r o)(x\RVB)} = 
0.971 [64] is very close to its maximal value, (s) = 1. We note that both the 
values of the overlap and of the average sign are much better than those ob- 
tained by a wave function that describes a magnetically ordered state, despite 
the smaller number of variational parameters |67j . Although the short-range 
RVB state is a very good variational ansatz, the bosonic representation of 
this state is rather difficult to handle in large clusters. Its systematic improve- 
ment by the inclusion of long-range valence bonds leads to a very severe sign 
problem, even at the variational level [55] . In this respect, following the rules 
discussed in section it is possible to obtain a fermionic representation of 
the short-range RVB state. The signs of the pairing function Jr^ri are given 
in Fig. ll.8l for open boundary conditions. Remarkably, this particular pattern 
leads to a 2 x 1 unit cell, which cannot be eliminated by using local SU(2) 
transformations of the type discussed in section 11.21 The variational RVB 
wave function is obtained by projecting the ground state of the BCS Hamil- 
tonian, with a particular choice of the couplings: the only nonzero parameters 
are the chemical potential \i and the nearest-neighbor singlet gap Ar^ri , in 
the limit — /x 3> |Ar,.R'| (so that the pairing function is proportional to the 
superconducting gap). The amplitude of the gap = A is uniform, 

while the appropriate phases are shown in Fig. 11.81 The BCS Hamiltonian is 
defined on a 2 x 1 unit cell and, therefore, is not translationally invariant. 
Despite the fact that it is invariant under an elementary translation T2 in the 
T2 = (1/2, y/3/2) direction, it is not invariant under an elementary translation 
T\ in the n = (1,0) direction. Nevertheless, this symmetry is recovered after 
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the projection Pq, making \pBCS) translationally invariant. Indeed, one can 
combine the translation operation Ti with the SU(2) gauge transformation 

for R — m\T\ + TO2T2 with 777,2 odd. Under the composite application of the 
transformations T\ and (|1.42[) . the projected BCS wave function does not 
change. Because the gauge transformation acts as an identity in the physical 
Hilbert space with singly occupied sites, \pBCS) is translationally invariant. 

Through this more convenient representation of the short-range RVB state 
by the pBCS wave function, it is possible to calculate various physical quan- 
tities using the standard variational Monte Carlo method. One example is 
the very accurate estimate of the variational energy per site in the thermo- 
dynamic limit, E/J = —0.5123(1) [64] ■ Another important advantage of the 
fermionic representation is that it is easy to improve the variational ansatz 
in a systematic way. The variational energy can be improved significantly by 
simply changing the chemical potential /1 from a large negative value to zero, 
see Table 11.21 We note that in this case \pBCS) is equivalent to a Gutzwiller- 
projected free fermion state with nearest-neighbor hoppings defined in a 2 x 1 
unit cell, because, through the SU(2) transformation of Eq (|1.8p . the off- 
diagonal pairing terms are transformed into kinetic terms. Further, the BCS 
Hamiltonian may be extended readily to include long-range valence bonds 
by the simple addition of nonzero Ar^ri or tR Ri terms. It is interesting to 
note that, within this approach, it is possible to obtain a variational energy 
E/J = —0.5357(1) lower than that obtained by starting from a magnetically 
ordered state and considered in Ref. [55] , see Table 11.21 

Finally, projected states have been also used to describe the ground state 
of the Hciscnberg Hamiltonian on the kagome lattice [69][70]. In this case, 
different possibilities for the mean-field Hamiltonian have been considered, 
with no BCS pairing but with non-trivial fluxes through the triangles and the 
hexagons of which the kagome structure is composed. In particular, the best 
variational state in this class can be found by taking 

H M F = - ^2 XR,R'C R ,aCR'^ + H.c, (1.43) 

{R,R').cr 



Table 1.2. Variational energy estimated in the thermodynamic limit for the anti- 
ferromagnetic Heisenberg model on the isotropic triangular lattice (J' = J). 



wave function 


E/J 


short-range RVB 


-0.5123(1) 


RVB with n = 


-0.5291(1) 


best RVB [Hi] 


-0.5357(1) 


BCS+Neel [68] 


-0.532(1) 
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with all the hoppings xr,R' having the same magnitude and producing a zero 
flux through the triangles and 7r flux through the hexagons. One may fix a 
particular gauge in which all Xrr 1 arc reai J see Fig. 11.91 In this gauge, the 
mean-field spectrum has Dirac nodes at k = (0,±7r/v3), and the variational 
state describes a U(l) Dirac spin liquid. Remarkably, this state should be sta- 
ble against dimcrization (i.e., it has a lower energy than simple valence-bond 
solids), in contrast to mean-field results [71] . Another competing mean-field 
state [7T] , which is obtained by giving the fermions chiral masses and is char- 
acterized by a broken time-reversal symmetry (with 8 flux through triangles 
and 7r — 8 flux through hexagons) , is also found to have a higher energy than 
the pure spin-liquid state. In this context, it would be valuable to compare 
the wave function proposed in Ref. [69] with the systematic improvement of 
the short-range RVB state which has a simple fermionic representation (see 
Fig.HU). 

1.8 Conclusions 

In summary, we have shown that projected wave functions containing both 
electronic pairing and magnetism provide an extremely powerful tool to study 
highly frustrated magnetic materials. In particular, these pBCS states may 
describe all known phases in one-dimensional systems, giving very accurate 
descriptions when compared to state-of-the-art DMRG calculations. Most im- 
portantly, variational wave functions may be easily generalized to treat higher 
dimensional systems: here we have presented in detail the case of the two- 
dimensional J\ — Ji model, as well as some examples of other frustrated lattices 
which have been considered in the recent past. 

The great advantage of this variational approach in comparison with other 
methods, such as DMRG, is that it can offer a transparent description of the 
ground-state wave function. Furthermore, the possibility of giving a physical 
interpretation of the unprojected BCS spectrum Ek, which is expected to be 
directly related to the true spin excitations, is very appealing. We demon- 
strated that this correspondence works very well in one dimension, both for 
gapless and for dimerized phases. In two dimensions, the situation is more 
complicated and we close by expressing the hope that future investigations 
may shed further light one the fascinating world of the low-energy properties 
of disordered magnetic systems. 
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